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Abstract. We investigate the phase diagram of the three-dimensional Hubbard model at half filling using 
quantum Monte Carlo (QMC) simulations. The antiferromagnetic Neel temperature Tjv is determined from 
the specific heat maximum in combination with finite-size scaling of the magnetic structure factor. Our 
results interpolate smoothly between the asymptotic solutions for weak and strong coupling, respectively, in 
contrast to previous QMC simulations. The location of the metal-insulator transition in the paramagnetic 
phase above Tjv is determined using the electronic compressibility as criterion. 
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1 Introduction 

The canonical lattice model for correlated electrons is the 
Hubbard model defined by the hamiltonian 



(1) 



<ij>,<r 



where cj^ (c itT ) are electron creation (annihilation) oper- 
ators, (ij) denotes a pair of neighboring lattice sites and 
t and U are the hopping matrix element and the onsite 
Coulomb energy, respectively. At half filling the ground 
state of the 3D Hubbard model on a simple cubic lat- 
tice has antiferromagnetic long range order for all posi- 

: tive values of U due to the perfect nesting of the Fermi 
surface. At finite temperatures the sublattice magnetiza- 
tion is reduced by thermal fluctuations and a transition 

, to a paramagnetic phase occurs at the Neel temperature 
Tn- In the strong coupling limit the energy scale for mag- 
netism is obtained by mapping the Hubbard model to 
the antiferromagnetic spin 1/2 Heisenberg model with ex- 
change coupling J = At 2 /U. For small U antifcrromag- 
netism results from the Fermi surface instability and the 
relevant energy scale is set by the BCS like expression 
Tm oc texp(— 1/{pqU)) where po is the density of states at 
the band center. 

There have been many attempts to calculate T/v over 
the whole range of U using both QMC simulations |^,|| 
and analytical approaches including variational methods 
Hartree Fock theory JtJ, strong coupling expansions 
j§, dynamical mean field theory (DMFT) the 
two-particle self-consistent formalism Jl2| and spin fluc- 
tuation theory |13|. A critical comparison of the various 



approaches can be found in refs. || and |12[ . A com- 
mon drawback of mean-field like approximations is that 
in the strong coupling limit they fail to reproduce the cor- 
rect Heisenberg result T/v = 3.83 t 2 /U 14J] and instead 
reduce to the Weiss molecular field theory with = 
6t 2 /U. Strong coupling expansions, on the other hand, 
break down in the Fermi surface instability regime. As to 
numerical methods, the results of previous QMC simu- 
lations do not agree with each other and have been 
questioned by Hasegawa |l5) who argued that they overes- 
timate Tjv considerably. Nevertheless, for lack of alterna- 
tives the QMC results of refs. have served as bench- 
marks for analytical approaches over the last ten years de- 
spite the controversy concerning their reliability. Clearly, 
there is need for improved QMC simulations of the 3D 
Hubbard model with better statistics and for larger sys- 
tems than previously accessible. 

Besides antiferromagnetism the second important phe- 
nomenon described by the half-filled Hubbard model is 
the interaction-induced metal-insulator transition, known 
as Mott-Hubbard transition fujfl . This transition occurs 
when the ratio of the interaction strength U and the band- 
width W exceeds some critical value of order one. Unfor- 
tunately, the presence of antiferromagnetic order makes it 
impossible to observe the Mott-Hubbard transition in the 
3D Hubbard model at temperatures below Tjy DMFT 
calculations for the fully frustrated Hubbard model (jl7) 
where antiferromagnetism is completely suppressed pre- 
dict a first order metal-insulator transition line that per- 
sists up to a critical point at (U c , T c ) in the phase diagram, 
followed by a crossover region above this critical point. In 
this paper we examine if such a first order transition line 
can also be observed in the 3D Hubbard model or if it is 
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Fig. 1. Energy as function of temperature for L — 4 and 
U — 4,6,8, 10, 12 (from bottom to top). Error bars are much 
smaller than the size of the symbols. The curves are obtained 
with the procedure described in the text. 
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Fig. 2. Specific heat as function of temperature for [7 = 6 and 
L — 4,6, 8, 10. The inset shows the position of the maximum 
of the specific heat vs. 1/L. 



completely occluded by the antiferromagnetic phase below 
T N . 



2 Magnetic phase transition 

We have studied the 3D Hubbard model on a simple cu- 
bic lattice with periodic boundary conditions using a fi- 
nite temperature, grand canonical QMC algorithm fll|[l9| 
which is based on a discrete Hubbard-Stratonovich decou- 
pling of the Hubbard interaction term. In this algorithm, 
the inverse temperature or imaginary time (3 has to be 
divided into a finite number of steps At which introduces 
an error oc AT 2 tU into the calculations. We have chosen 
At 2 OJ — 0.1 after making sure that the results are not 
significantly affected by the extrapolation At — > 0. A typ- 
ical simulation on a L x L x L lattice consisted of between 
100000 (L = 4) and 2000 (L = 10) measurement sweeps 
that were grouped into 20 blocks in order to estimate the 
statistical error of the QMC data. QMC simulations fre- 
quently suffer from the minus sign problem, i.e. the fact 
that the fermionic determinant that serves as probabil- 
ity weight function is not always positive, in particular at 
low temperatures. At half filling where most of our sim- 
ulations were performed this problem does not exist due 
to particle-hole symmetry, but even in the calculations 
of the compressibility that require simulations away from 
half filling we never encountered any serious minus sign 
problem since we worked at relatively high temperatures. 

Our first goal is to determine the antiferromagnetic 
Neel temperature as function of the interaction strength 
U. To this end we have calculated the specific heat which 
is the central quantity to characterize thermodynamical 
properties of many-particle systems and is easily accessi- 
ble both experimentally and from numerical simulations. 
Since phase transitions are usually accompanied by a max- 
imum or a divergence of the specific heat the transition 
temperature can be determined from specific heat data 



without measuring the order parameter directly. Of course, 
not every maximum of the specific heat is associated with 
a phase transition; one has to verify that the order param- 
eter indeed changes from zero to a nonzero value at the 
presumed transition temperature. In the Hubbard model 
for sufficiently large U there are two maxima of the spe- 
cific heat, one at low temperatures cat 2 /U reflecting spin 
excitations, the second on a higher energy scale oc U as- 
sociated with charge excitations. Here we concentrate on 
the low temperature peak which we use to determine the 
Neel temperature. 

The specific heat can either be calculated via C = 
dE/dT or from the energy fluctuations, C oc (H 2 ) — (H) 2 . 
We employ the first method which turns out to be numeri- 
cally more accurate. The following procedure is used: First 
the energy is calculated with QMC in a temperature range 
where the phase transition is expected to take place. The 
simulations are performed at fixed imaginary time slice 
At for all temperatures. An extrapolation A — > is not 
necessary since the finite At corrections to the energy are 
nearly temperature independent over the relatively small 
temperature range under consideration and even taking 
into account corrections to linear order in T does not affect 
the position of the specific heat maximum. In order to cal- 
culate the specific heat from the energy data the tempera- 
ture interval is divided in a number of equidistant values. 
For each of these temperatures Tj the value of the specific 
heat Ci is determined such, that i) the resulting curve is 
as smooth as possible and ii) the energies E^ obtained 
by numerical integration of the specific heat values Ci are 
as close as possible to the QMC data. In practice this is 
achieved by minimizing the quantity XiA + X2B where 
A = 5D 4 (C<+i — 2Ci + Ci-1) 2 controls the smoothness of 
the fit while B = £((£fit - ^QMCV^QMC) 2 measures 
the deviation between the fit and the QMC data. This pro- 
cedure avoids to chose a specific functional form for the fit 
function which could bias the results in some way or the 
other. The ratio of the parameters Ai and A2 is fixed such 
that the energy values obtained by the minimization pro- 
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cedure lie within the statistical errors <JQy[(j of the QMC 
data. We have checked that the final results are very ro- 
bust against moderate variations of the ratio A1/A2. In 
order to estimate the error bars for the specific heat we 
have performed an average over many different input en- 
ergy data sets obtained by adding random noise of the or- 



der of a 



QMC 



to the mean values E, 



QMC 



Of course, any 



noise-reducing algorithm tends to wash out sharp struc- 
tures like cusps or discontinuities. On the other hand, the 
specific heat curves obtained from extensive QMC simula- 
tions of the three-dimensional spin 1/2 Heisenberg model 
p?| are as smooth as ours indicating that there exist no 
sharp structures that could be lost by the procedure that 
we use to calculate the specific heat, at least for the sys- 
tem sizes that we consider. In any case, the location of the 
maximum which we are mainly interested in is only very 
little affected. 

Fig. [l] shows the energy as function of temperature 
for L = 4 and several values of U. Here and in the fol- 
lowing figures we have fixed the energy scale by setting 
t = 1. The curves connecting the data are obtained with 
the fitting procedure described above. They all have a 
point of inflection indicating that the specific heat max- 
ima are contained in the respective temperature intervals. 
In Fig. |^ the specific heat for U = 6 is displayed as func- 
tion of temperature for L = 4, 6, 8 and 10. The maximum 
value of C increases somewhat with increasing system size 
while simultaneously the position of the maximum T max is 
slightly shifted to lower temperatures. There is no indica- 
tion of a divergence for large systems. This is in agreement 
with the assumption that the half-filled Hubbard model 
belongs to the 3D Heisenberg universality class with a 
negative exponent a ~ —0.1 1 which means that the spe- 
cific heat has only a cusp but no divergence at the criti- 
cal temperature. High precisio n Q MC simulations of the 
3D spin 1/2 Heisenberg model EQ] confirm this behavior. 
In the Heisenberg model the shift of the maximum of C 
between L = 4 and the infinite lattice is about five per- 
cent. In the inset of Fig. || the peak temperature T max 
is plotted vs 1/L indicating that in the Hubbard model 
finite-size corrections are quite small as well. A linear fit 
yields T N = 0.31 ± 0.01 . 

We now demonstrate that the peaks in the specific heat 
are indeed associated with the antiferromagnetic phase 
transition. To this end we have calculated the magnetic 
structure factor 

= J3 E e lQ(R '- R ^ (fat ~ ^)(%T ~ W) ( 2 ) 

where Q = (71", it, n) is the antiferromagnetic wave vector. 
S(Q) is related to the sublattice magnetization m via 



S(Q) 

T 3 



f(L) 



(3) 



where f(L) — > for L — > 00. In order to extrapolate the 
finite lattice QMC data to the thermodynamic limit one 
has to know the asymptotic behavior of f(L) for large 
L. Right at the critical temperature the structure factor 




Fig. 3. Finite size scaling of the magnetic structure factor for 
(7 = 6 and T = 0.3, 0.36 and 0.4 (from top to bottom). 



S(Q) should scale with the system size as L 2 ~ v where 
r\ k, 0.03 for the 3D Heisenberg universality class p(J and 
therefore f(L) cx L~ 103 . In the ordered phase at low tem- 
peratures, spin wave theory predicts /(T) cx L~ 2 assum- 
ing a linear magnon dispersion. On the other hand, in the 
paramagnetic phase above TV spin correlations decay ex- 
ponentially and we expect f(L) cx L~ 3 provided the cor- 
relation length is smaller than the system size. In order to 
take into account all these possibilities we have extrapo- 
lated our QMC data using f(L) cx where A itself is 
a fit parameter. Since the asymptotic behavior of f(L) is 
only reached for sufficiently large lattices we have checked 
that the omission of the data for small lattices (L = 4 and 
L = 6) does not lead to different conclusions concerning 
the existence or absence of antiferromagnetic long range 
order. 

The results of such a finite-size extrapolation are dis- 
played in Fig. I for U = 6 and T = 0.3, 0.36 and 0.4. While 
the curves for the two higher temperatures extrapolate to 
a value close to zero indicating the absence of long-range 
order, the curve for T = 0.3 yields a finite value corre- 
sponding to a finite sublattice magnetization. This behav- 
ior is in agreement with the value TV ps 0.31 that we have 
extracted from the specific heat data. Preliminary results 
concerning the magnetic structure factor for other values 
of U have been published elsewhere J2lj . A more compre- 
hensive presentation including computational details can 
be found in ||. 

We summarize our results concerning the magnetic 
phase transition in Fig. || where the Nccl temperature of 
the 3D Hubbard model obtained from different methods 
is displayed as function of U. For comparison the asymp- 
totic behavior in the weak and strong coupling limit is 
also shown. In Hartree Fock theory, applicable for small 
U, Tjv is determined by the gap equation 



2 f , p(e) , e 
U = J e~ * 2TV 



(4) 



where p(e) is the density of states. As pointed out by van 
Dongen W the true asymptotic Neel temperature in the 
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Fig. 4. Magnetic phase diagram of the half-filled Hubbard 
model. Neel temperature Tn as function of U from various ap- 
proaches: QMC, this work (filled circles; dots are meant as 
a guide to the eye only), QMC Q (open triangles), QMC 
H (filled triangles), DMFT § (open circles), DMFT @ 
(open squares), modified Hartree-Fock theory |7f (solid curve), 
Heisenber g li mit from high temperature expansions, Tn = 
3.83t 2 /(7 114] (dashed curve), Weiss molecular field theory, 
T™ f = 6t 2 /U (dotted curve). 



weak coupling limit is reduced by a factor q ss 0.282 com- 
pared to the solution of (|j). We have included this re- 
duction factor in the curve shown in the figure. In the 
strong coupling limit the Hubbard model can be mapped 
to the spin 1/2 Heisenberg model where the critical tem- 
perature is known from high temperature series |l4| and 
QMC simulations @ which yield T N = 3.83t 2 /U, com- 
pared to the Weiss molecular field result Tjy = 6t 2 /U. 
Our QMC results interpolate smoothly between weak and 
strong coupling asymptotics whereas the old QMC data 
of refs. [|J and Q are clearly off in both limits. 

Comparing our results with DMFT data is not straight- 
forward since in DMFT calculations mostly a Gaussian or 
a semi-elliptic density of states is used which both dif- 
fer from the non-interacting density of states of the 3D 
Hubbard model. To convert energies we have expressed 
the hopping matrix element t (which is our energy unit) 
in terms of the second moment of the density of states 
via t = \J (e 2 )/6. In practice this means that the energy 
data from ref. ||] have been multiplied by a factor 2\/3 
and those from ref. [[ll| by a factor of 2\f%. In the weak 
coupling regime U < 6 there is good agreement between 
the DMFT results and our QMC data while for intermedi- 
ate and large values of U the DMFT yields substantially 
higher values of T/y. This is not surprising since in the 
limit U — + oo the DMFT reduces to the Weiss molecular 
field theory of the Heisenberg model. There are however 
significant discrepancies between the DMFT data from 
ref. p and ref. [jll] for large U. This might be due to the 
fact that in the first case a Gaussian and in the latter a 
semi-elliptic density of states was used. It should also be 
noted that the DMFT data from approach the molec- 
ular field limit from above while the curve of 10] and also 
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Fig. 5. Particle density n as function of the chemical potential 
H for U = 3, 4, 5, 6, 7, 8, 9, 10, 12 (from top to bottom), T = 0.45 
and L = 4. The curves are polynomial fits. 

our QMC data lie always below their respective strong 
coupling asymptote. 



3 Metal insulator transition 

At a true metal insulator transition the DC conductivity 
a drops to zero when some control parameter is changed 
across a critical value. This is strictly speaking only possi- 
ble at zero temperature whereas for T > the conductiv- 
ity remains always finite due to thermal activation. Un- 
fortunately, the DC conductivity - although quite easily 
accessible in experiment - is very hard to obtain from 
QMC simulations, since the required analytic continuation 
of the current-current correlation function from Matsub- 
ara to real frequencies is numerically a very hard prob- 
lem when the input data are noisy. We have therefore 
based our considerations on calculations of the electronic 
compressibility k = dnjd[i which can be obtained with 
high accuracy from QMC simulations. Although there is 
no simple relation between a and k both are expected to 
be finite in the metallic and exponentially small in the 
insulating regime. 

DMFT calculations employing the iterated perturba- 
tion theory Jl7| have revealed the following scenario for the 
Mott Hubbard transition. In the fully frustrated Hubbard 
model there exists a first-order metal-insulator transition 
line due to a coexistence regime of metallic and insulating 
solutions at low temperatures. This first order line ends 
at a critical point, reminiscent of an ordinary liquid gas 
transition |23| . Recent numerical work has corroborated 
this scenario. In fact, the first order transition can also be 
observed as a jump in the average number of doubly occu- 
pied sites |^3| and it is expected that the compressibility 
is discontinuous as well. 

Fig. shows how the particle density n depends on the 
chemical potential fi for T — 0.45 which is slightly above 
the antiferromagnetic phase boundary. Due to particle- 
hole symmetry n— 1 is an odd function of fi—U/2. We have 
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Fig. 6. Electronic compressibility n as function of U. The curve 
is a polynomial fit for the range 3 < U < 8. 



therefore used the polynomial f(x) = ax + bx 3 



to fit 



the QMC data. The coefficient a yields the compressibility 
k which is displayed as function of U in Fig. ^|. There is no 
sign of a discontinuity as should be observed in the case 
of a first order transition. Instead the compressibility de- 
creases smoothly up to U ~ 8 and can be very accurately 
approximated by a second order polynomial in this range, 
as indicated by the dashed curve in the figure. Afterwards 
it turns over into an exponential-like tail that we associate 
with the crossover regime where only thermal activation 
across the gap contributes to the compressibility. Unfor- 
tunately our data are not precise enough to extract the 
U dependence of the gap. The observation of a crossover 
regime instead of a true metal-insulator transition is in 
agreement with the location of the critical point calcu- 
lated within DMFT. For a semi-elliptic density of states 
the critical values are U c ~ 11.7 and T c ~ 0.127 |23|, con- 
verted to our units which is by a factor of more than two 
below the maximum Neel temperature i.e. deep inside the 
antifcrromagnetic phase. 



4 Conclusions 

We have performed QMC simulations of the half-filled 3D 
Hubbard model on simple cubic lattices of up to 10 3 sites. 
Using finite-size scaling of the magnetic structure factor 
it is shown that the low-temperature maximum of the 
specific heat coincides with the antiferromagnetic phase 
transition. The Neel temperature thus obtained interpo- 
lates smoothly between the analytic solutions for weak 
and strong coupling but differs significantly from the re- 
sults of previous QMC simulations. The shape of the spe- 
cific heat curves close to the phase transition is similar to 
the one obtained for the 3D spin 1/2 Heisenberg model 
using high precision QMC simulations |^0| indicating that 
both models are in the same universality class. There is 
no indication of a first-order transition for small values of 
U contrary to a conjecture established in previous QMC 
simulations 0,0]. 



In order to investigate the transition from metallic 
to insulating behavior in the paramagnetic phase above 
Tjv we have calculated the electronic compressibility. The 
transition appears to be broadened in accordance with 
the crossover scenario developed in the framework of the 
DMFT. For T — 0.45 the crossover regime extends over 
the range 9 < U < 11 which is somewhat below what is 
obtained in DMFT calculations (||. 

It would be very interesting to perform QMC simula- 
tions for a Hubbard model where antiferromagnetism is 
suppressed by adding a next-nearest neighbor hopping t' 
in order to confirm the existence of the first order tran- 
sition line obtained within the DMFT. It is however to 
be feared that in this case the notorious minus sign prob- 
lem will prohibit efficient QMC simulations at low enough 
temperatures. We estimate that the average sign of the 
fermion determinant behaves as cx exp(— f3t') which means 
that for temperatures of the order of t' the minus sign 
problem becomes severe and therefore the development of 
improved QMC algorithms is needed to study this prob- 
lem. 

We thank P. van Dongen and P. Schwab for useful 
discussions. This work was supported by the DFG as a 
part of the project metal-insulator transition and mag- 
netism in highly correlated transition-metal chalcogenids 
(HO 955/2). 
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